// Copyright (c) Lawrence Livermore National Security, LLC and other VisIt
// Project developers.  See the top-level LICENSE file for dates and other
// details.  No copyright assignment is required to contribute to VisIt.

// ************************************************************************* //
//                              avtTopologyFilter.C                          //
// ************************************************************************* //

#include <avtTopologyFilter.h>

#include <vtkCell.h>
#include <vtkCellData.h>
#include <vtkCellDataToPointData.h>
#include <vtkContourFilter.h>
#include <vtkDataSet.h>
#include <vtkExecutive.h>
#include <vtkFloatArray.h>
#include <vtkIdList.h>
#include <vtkPointData.h>
#include <vtkPolyData.h>
#include <vtkRectilinearGrid.h>
#include <vtkStructuredGrid.h>

#include <ImproperUseException.h>

#include <math.h>

#include <string>
#include <vector>

using std::string;
using std::vector;

const string avtTopologyFilter::labelNames[4] =
    { "min", "min plateau", "max plateau", "max" };


// ****************************************************************************
//  Method: avtTopologyFilter constructor
//
//  Programmer: haddox1 -- generated by xml2info
//  Creation:   Thu Jun 26 15:59:09 PST 2003
//
// ****************************************************************************

avtTopologyFilter::avtTopologyFilter()
{
    tolerance = 1e-6;
    hitpercent = 0.6;
    
    cf = vtkContourFilter::New();
    cf->SetUseScalarTree(1);
    cf->SetNumberOfContours(1);
}


// ****************************************************************************
//  Method: avtTopologyFilter destructor
//
//  Programmer: haddox1 -- generated by xml2info
//  Creation:   Thu Jun 26 15:59:09 PST 2003
//
// ****************************************************************************

avtTopologyFilter::~avtTopologyFilter()
{
    if (cf)
        cf->Delete();
}


// ****************************************************************************
//  Method: avtTopologyFilter::ExecuteDataTree
//
//  Purpose:
//      Does the actual VTK code to modify the dataset. Returns one
//      dataset for each topology value.
//
//  Arguments:
//      in_dr     The input data representation.
//
//  Returns:      The output data representation.
//
//  Programmer: haddox1 -- generated by xml2info
//  Creation:   Thu Jun 26 15:59:09 PST 2003
//
//  Modifications:
//    Kathleen Bonnell, Tue May 16 13:50:50 PDT 2006
//    Removed call to ds->SetSource(NULL): with new vtk pipeline it also removes
//    necessary information form the dataset.  Changed SetOutput to 
//    GetExecutive()->SetOutputData.
//
//    Eric Brugger, Tue Aug 19 11:41:51 PDT 2014
//    Modified the class to work with avtDataRepresentation.
//
// ****************************************************************************

avtDataTree_p
avtTopologyFilter::ExecuteDataTree(avtDataRepresentation *in_dr)
{
    //
    // Get the VTK data set and the domain number.
    //
    vtkDataSet *in_ds = in_dr->GetDataVTK();
    int domain = in_dr->GetDomain();

    int jump;
    
    if (in_ds->GetDataObjectType() == VTK_RECTILINEAR_GRID)
        jump = ((vtkRectilinearGrid*)(in_ds))->GetDimensions()[0];
    else
        jump = -1;

    int npts   = in_ds->GetNumberOfPoints();

    //
    // Get (or make) the array of node based scalar data.
    //
    bool ownDataArray = false;
    vtkDataArray *data = in_ds->GetPointData()->GetScalars();

    if (!data)
    {
        ownDataArray = true;
        if (!in_ds->GetCellData()->GetScalars())
            EXCEPTION1(ImproperUseException, "only for scalar data");
        
        vtkCellDataToPointData *myFilter = vtkCellDataToPointData::New();
        myFilter->SetInputData(in_ds);
        myFilter->Update();
        data = myFilter->GetOutput()->GetPointData()->GetScalars();
        data->Register(NULL);
        myFilter->Delete();
    }
    
    //
    // Allocate the arrays for storing the four types of values to contour.
    //
    vtkFloatArray *rv[4];
    int i;
    for (i = 0; i < 4; ++i)
    {
        rv[i] = vtkFloatArray::New();
        rv[i]->SetNumberOfTuples(npts);
    }
    
    bool is3D = in_ds->GetBounds()[4] != in_ds->GetBounds()[5];
    
    if (is3D)
        EXCEPTION1(ImproperUseException, "Currently only for 2D datasets.");
    
    //
    // If we're a rectilinear grid, we can calculate node based gradients,
    // and work directly off of those.
    //
    vtkDataArray *gradients = NULL;
    if (jump != -1)
       gradients = GetGradient(in_ds, data);
    
    for (i = 0 ; i < npts ; i++)
    {
        bool sideValid = true;
        bool topValid = true;

        double ptCoords[3];
        in_ds->GetPoint(i, ptCoords);

        //
        // Get the gradients around this point. Also find if the sides or
        // vertical points are invalid.
        //
        double up, down, right, left;
        double leftBlock[3], rightBlock[3], upBlock[3], downBlock[3];

        if (jump == -1)
        {
            sideValid = GetLocalGradient(0, in_ds, data, leftBlock, i);
            if (GetLocalGradient(1, in_ds, data, rightBlock, i) == false)
                sideValid = false;;
            topValid = GetLocalGradient(2, in_ds, data, upBlock, i);
            if (GetLocalGradient(3, in_ds, data, downBlock, i) == false)
                topValid = false;
        }
        else
        {
            gradients->GetTuple((i - 1 + npts) % npts, leftBlock);
            gradients->GetTuple((i + 1) % npts, rightBlock);
            gradients->GetTuple((i + jump) % npts, upBlock);
            gradients->GetTuple((i - jump + npts) % npts, downBlock);

            if (i < jump || i >= npts - jump)
               topValid = false;
            if (i % jump == 0 || (i + 1) % jump == 0)
               sideValid = false;
        }
        
        left = leftBlock[0];
        right = rightBlock[0];
        up = upBlock[1];
        down = downBlock[1];

        //
        // Convenience booleans:
        //   t*         :   Gradient is coming towards this node
        //   a*         :   Gradient is going away from this node
        //   z*         :   Gradient is zero (based on tolerance)
        //   [vh]S      :   Gradient is significant (based on hitpercent)
        //                  for two gradients [vertical and horizontal]
        //   [udrl]S    :   Gradient is significant (based on hitpercent)
        //                  for one gradient.
        //
        
        bool tu = up < 0;
        bool td = down > 0;
        bool tr = right < 0;
        bool tl = left > 0;
        bool au = up > 0;
        bool ad = down < 0;
        bool ar = right > 0;
        bool al = left < 0;

        double validThreshold = 2 * hitpercent;

        double upPercent = fabs(upBlock[1]) / 
                            (fabs(upBlock[0]) + fabs(upBlock[1]));
        double downPercent = fabs(downBlock[1]) / 
                            (fabs(downBlock[0]) + fabs(downBlock[1]));
        double leftPercent = fabs(leftBlock[0]) / 
                            (fabs(leftBlock[0]) + fabs(leftBlock[1]));
        double rightPercent = fabs(rightBlock[0]) / 
                            (fabs(rightBlock[0]) + fabs(rightBlock[1]));

        bool vS = upPercent + downPercent >= validThreshold;
        bool hS = leftPercent + rightPercent >= validThreshold;

        bool uS = upPercent >= hitpercent;
        bool dS = downPercent >= hitpercent;
        bool rS = rightPercent >= hitpercent;
        bool lS = leftPercent >= hitpercent;

        double epsilon = tolerance;
        bool zu = fabs(upBlock[0]) < epsilon && fabs(upBlock[1]) < epsilon;
        bool zd = fabs(downBlock[0]) < epsilon && fabs(downBlock[1]) < epsilon;
        bool zl = fabs(leftBlock[0]) < epsilon && fabs(leftBlock[1]) < epsilon;
        bool zr = fabs(rightBlock[0]) < epsilon && fabs(rightBlock[1]) < epsilon;
        
        // Default all the values to 0.
        rv[0]->SetTuple1(i, 0);
        rv[1]->SetTuple1(i, 0);
        rv[2]->SetTuple1(i, 0);
        rv[3]->SetTuple1(i, 0);

        if (is3D)
        {            
            ;
        }
        else
        {
            //
            // Check for saddle points to disqualify them
            //
            if ((tu && td && ar && al && vS && hS)
             || (au && ad && tr && tl && vS && hS))
                ;
            
            // Check for constant area first
            else if (zu && zd && zr && zl)
                ;

            // Find gradient to zero
            else if (tu && zd && uS && topValid)
                rv[2]->SetTuple1(i, 1);
            else if (td && zu && dS && topValid)
                rv[2]->SetTuple1(i, 1);
            else if (tl && zr && lS && sideValid)
                rv[2]->SetTuple1(i, 1);
            else if (tr && zl && rS && sideValid)
                rv[2]->SetTuple1(i, 1);

            else if (au && zd && uS && topValid)
                rv[1]->SetTuple1(i, 1);
            else if (ad && zu && dS && topValid)
                rv[1]->SetTuple1(i, 1);
            else if (al && zr && lS && sideValid)
                rv[1]->SetTuple1(i, 1);
            else if (ar && zl && rS && sideValid)
                rv[1]->SetTuple1(i, 1);
            
            // Find max ridges
            //  Max ridges are where a pair of oppossing vectors have that
            //  opposite sign.

            // Check top / bottom:
            //  This says that the y components of grads above and below
            //  are pointing towards me, and that the y component is
            //  significant.
            else if (tu && td && vS && topValid)
                rv[3]->SetTuple1(i, 1);
            else if (tr && tl && hS && sideValid)
                rv[3]->SetTuple1(i, 1);

            // Find min ridges:
            else if (au && ad && vS && topValid)
                rv[0]->SetTuple1(i, 1);
            else if (ar && al && hS && sideValid)
                rv[0]->SetTuple1(i, 1);
        }
    }
    
    //
    // Now we contour the four variables to draw the lines.
    //
    
    vtkDataSet *contourInput = (vtkDataSet *)in_ds->NewInstance();
    contourInput->ShallowCopy(in_ds);

    // Make this variable? 
    float contourLevel = 0.999;
    
    vector<vtkPolyData *> contours;
    vector<string> labels;
    
    size_t cl;
    for (cl = 0; cl < 4; ++cl)
    {
        contourInput->GetPointData()->SetScalars(rv[cl]);
        cf->SetInputData(contourInput);
        cf->SetValue(0, contourLevel);
        cf->Update();
        vtkPolyData *pd = vtkPolyData::New();
        pd->ShallowCopy(cf->GetOutput());
        contours.push_back(pd);
        labels.push_back(labelNames[cl]);
    }

    cf->SetInputData(NULL);

    rv[0]->Delete();
    rv[1]->Delete();
    rv[2]->Delete();
    rv[3]->Delete();
    contourInput->Delete();
    if (gradients)
        gradients->Delete();
    if (ownDataArray)
        data->Delete();
    

    size_t ccount = contours.size();
    if (ccount == 0)
        return NULL;

    vtkDataSet **out_ds = new vtkDataSet*[ccount];
    for (cl = 0; cl < ccount; ++cl)
        out_ds[cl] = contours[cl];

    avtDataTree_p returnDT = new avtDataTree((int)ccount, out_ds, domain, labels);
    
    return returnDT;
}


// ****************************************************************************
//  Method: avtTopologyFilter::GetGradient
//
//  Purpose:
//      Get gradient of a scalar variable from a dataset centered on nodes.
//
//  Arguments:
//      in_ds               The input dataset.
//      scalarValues        The scalar variable (node based)
//
//  Returns:      The derived variable.  The calling class must free this
//                memory.
//
//  Programmer:   Akira Haddox 
//  Creation:     June 26, 2003
//
// ****************************************************************************

vtkDataArray *
avtTopologyFilter::GetGradient(vtkDataSet *in_ds, vtkDataArray *scalarValues)
{
    int nPoints = in_ds->GetNumberOfPoints();

    vtkDataArray *results = vtkFloatArray::New();
    results->SetNumberOfComponents(3);
    results->SetNumberOfTuples(nPoints);        
    
    for (int nodeId = 0 ; nodeId < nPoints; nodeId++)
    {
        float xDELTA=1e6, yDELTA=1e6, zDELTA=1e6;
        
        double nodeCoords[3];
        in_ds->GetPoint(nodeId, nodeCoords);
        float  nodeValue  = scalarValues->GetComponent(nodeId, 0);
        
        vtkIdList *neighborCellIds = vtkIdList::New();
        vtkIdList *myNodeId = vtkIdList::New();

        myNodeId->SetNumberOfIds(1);
        myNodeId->SetId(0, nodeId);
        
        in_ds->GetCellNeighbors(-1, myNodeId, neighborCellIds);

        myNodeId->Delete();

        int nCells=neighborCellIds->GetNumberOfIds();

        // Find appropriate deltas
        for (int ci = 0 ; ci < nCells ; ci++)
        {
            double *bounds = in_ds->GetCell(neighborCellIds->GetId(ci))
                                                                 ->GetBounds();
            if (bounds[1]-bounds[0] < xDELTA*5)
            {
                xDELTA = (bounds[1]-bounds[0]) / 5.0;
            }
            if (bounds[3]-bounds[2] < yDELTA*5)
            {
                yDELTA = (bounds[3]-bounds[2]) / 5.0;
            }
            if (bounds[5]-bounds[4] < zDELTA*5)
            {
                zDELTA = (bounds[5]-bounds[4]) / 5.0;
            }
        }
        
        double xComponent, yComponent, zComponent;
        
        xComponent=EvaluateComponent(nodeCoords[0],nodeCoords[1],nodeCoords[2],
                xDELTA, 0, 0, nodeValue, in_ds, scalarValues, neighborCellIds);
        yComponent=EvaluateComponent(nodeCoords[0],nodeCoords[1],nodeCoords[2],
                0, yDELTA, 0, nodeValue, in_ds, scalarValues, neighborCellIds);
        if (zDELTA == 0)
        {
            zComponent = 0;
        }
        else
        {
            zComponent = EvaluateComponent(nodeCoords[0],nodeCoords[1],
                                           nodeCoords[2], 0, 0, zDELTA,
                                           nodeValue, in_ds, scalarValues,
                                           neighborCellIds);
        }

        results->SetTuple3(nodeId, (float)xComponent, (float)yComponent,
                                   (float)zComponent);
    }
    
    return results;
}


// ****************************************************************************
//  Method: avtTopologyFilter::EvaluateComponent
//
//  Purpose:
//      Calculate the directional derivative from the given point in the given
//      deltas.
//
//  Arguments:
//    x,y,z           The point.
//    dx, dy, dz      The direction, in delta scale: It is indended for (2/3)
//                    of them to be 0.
//    value           The value at the point.
//    in_ds           The dataset.
//    scalarValues    The values for the dataset.
//    neighborCells   The cells to look in.
//
//  Returns:      The value at that point.
//
//  Programmer:   Akira Haddox 
//  Creation:     June 26, 2003
//
// ****************************************************************************

float avtTopologyFilter::EvaluateComponent(float x, float y, float z,
            float dx, float dy, float dz, float value, vtkDataSet *in_ds,
            vtkDataArray *scalarValues, vtkIdList *neighborCells)
{
    int deltaMultiplier = 2;
    float upper;
    bool  success = true;

    upper = EvaluateValue(x+dx, y+dy, z+dz, in_ds, scalarValues, neighborCells,
                          success);
    if (!success)
    {
        upper=value;
        --deltaMultiplier;
    }

    float lower;
    success = true;
    lower = EvaluateValue(x-dx, y-dy, z-dz, in_ds, scalarValues, neighborCells,
                          success);

    if (!success)
    {
        lower=value;
        --deltaMultiplier;
    }

    if (!deltaMultiplier)
    {
        return 0;
    }
    
    return (upper-lower)/ (double(deltaMultiplier)*(dx+dy+dz));
}


// ****************************************************************************
//  Method: avtTopologyFilter::EvaluateValue
//
//  Purpose:
//      Calculate the value at a point.
//
//  Arguments:
//    x,y,z           The point.
//    in_ds           The dataset.
//    scalarValues    The values for the dataset.
//    neighborCells   The cells to look in.
//    success         Set to false if there is a problem.
//
//  Returns:      The value at that point.
//
//  Programmer:   Akira Haddox 
//  Creation:     June 26, 2003
//
// ****************************************************************************

float avtTopologyFilter::EvaluateValue(float x, float y, float z, 
                               vtkDataSet *in_ds, vtkDataArray *scalarValues,
                               vtkIdList *neighborCells, bool &success)
{
    // Find which cell contains this point
    double coords[3] = {x,y,z};
    int   junk2;
    double junk3[3];
    double junk4;
    double weights[8];  // This needs to be the max number of points a cell has

    double *abnormalWeights = NULL; // In case of more than 8 points
    
    int cellId;
    vtkCell *c = NULL;
    for (cellId = 0 ; cellId < neighborCells->GetNumberOfIds() ; cellId++)
    {
        c = in_ds->GetCell(neighborCells->GetId(cellId));
    
        if (c->GetNumberOfPoints() > 8)
        {
            abnormalWeights = new double[c->GetNumberOfPoints()];
            if (c->EvaluatePosition(coords, NULL, junk2, junk3, junk4, 
                                    abnormalWeights) == 1)
            {
                break;
            }
            
            delete[] abnormalWeights;
            abnormalWeights=NULL;
        }
        
        else if (c->EvaluatePosition(coords, NULL, junk2, junk3, junk4, 
                                     weights) == 1)
        {
            break;
        }
    }

    if (cellId == neighborCells->GetNumberOfIds())
    {
        success = false;
        return 0.;
    }
    success = true;
    
    double value = 0.;
    if (abnormalWeights)
    {
        for (int k=0 ; k < c->GetNumberOfPoints() ; k++)
        {
            int pt = c->GetPointId(k);
            value += abnormalWeights[k] * scalarValues->GetComponent(pt,0);
        }
        delete [] abnormalWeights;
    }
    else
    {
        for (int k = 0 ; k < c->GetNumberOfPoints() ; k++)
        {
            int pt = c->GetPointId(k);
            value += weights[k] * scalarValues->GetComponent(pt, 0);
        }
    }
    return value;    
} 


// ****************************************************************************
//  Method: avtTopologyFilter::GetLocalGradient
//
//  Purpose:
//      Get gradient of a scalar variable from a dataset in a small step
//      away from a point.
//
//  Arguments:
//      which               left = 0, right = 1, up = 2, down = 3
//      in_ds               The input dataset.
//      data                The scalar variable (node based)
//      block               To return the gradient
//      int                 The relative node id 
//
//  Returns:      True on success. False on failure (commonly due to
//                being on the edge of a domain).
//
//  Programmer:   Akira Haddox 
//  Creation:     June 26, 2003
//
// ****************************************************************************

bool
avtTopologyFilter::GetLocalGradient(int which, vtkDataSet *in_ds,
                                    vtkDataArray *data,
                                    double block[3], int nodeId)
{
    vtkIdList *neighborCellIds = vtkIdList::New();
    vtkIdList *myNodeId = vtkIdList::New();

    myNodeId->SetNumberOfIds(1);
    myNodeId->SetId(0, nodeId);
        
    in_ds->GetCellNeighbors(-1, myNodeId, neighborCellIds);

    myNodeId->Delete();
    int nCells=neighborCellIds->GetNumberOfIds();
    
    double xDELTA=1e6, yDELTA=1e6, zDELTA=1e6;
    
    // Find appropriate deltas
    for (int ci = 0 ; ci < nCells ; ci++)
    {
        double *bounds = in_ds->GetCell(neighborCellIds->GetId(ci))
                                                             ->GetBounds();
        if (bounds[1]-bounds[0] < xDELTA*5)
        {
            xDELTA = (bounds[1]-bounds[0]) / 5.0;
        }
        if (bounds[3]-bounds[2] < yDELTA*5)
        {
            yDELTA = (bounds[3]-bounds[2]) / 5.0;
        }
        if (bounds[5]-bounds[4] < zDELTA*5)
        {
            zDELTA = (bounds[5]-bounds[4]) / 5.0;
        }
    }

    //
    // An iterative approach. If our point fails to land in one of the
    // neighboring cells, we pull back our step amount and try again.
    // After a few tries (and sharp pulling back), we give up, concluding
    // that the point is probably outside of the dataset.
    //
    const int attempts = 4;
    int run;
    for (run = 0; run < attempts; ++run)
    {
 
        double nc[3];
        in_ds->GetPoint(nodeId, nc);

        if (which == 0)
            nc[0] -= xDELTA * 2;
        if (which == 1)
            nc[0] += xDELTA * 2;
        if (which == 2)
            nc[1] += yDELTA * 2;
        if (which == 3)
            nc[1] -= yDELTA * 2;

        bool succ = true;
        float val = EvaluateValue(nc[0], nc[1], nc[2], in_ds, data, neighborCellIds, succ);

        if (!succ)
        {
            xDELTA /= 3.0;
            yDELTA /= 3.0;
            zDELTA /= 3.0;
            continue;
        }
    
        block[0] = EvaluateComponent(nc[0], nc[1], nc[2], xDELTA, 0, 0, val,
                                     in_ds, data, neighborCellIds);
        block[1] = EvaluateComponent(nc[0], nc[1], nc[2], 0, yDELTA, 0, val,
                                     in_ds, data, neighborCellIds);
        block[2] = 0;
        return true;
    }

    // We failed. Must be outside the dataset.
    block[0] = block[1] = block[2] = 0;
    return false;
}


// ****************************************************************************
//  Method: avtTopologyFilter::UpdateDataObjectInfo
//
//  Purpose:
//      Allows the filter to change its output's data object information, which
//      is a description of the data object.
//
//  Programmer: haddox1 -- generated by xml2info
//  Creation:   Thu Jun 26 15:59:09 PST 2003
//
// ****************************************************************************

void
avtTopologyFilter::UpdateDataObjectInfo(void)
{
    GetOutput()->GetInfo().GetValidity().InvalidateZones();
    GetOutput()->GetInfo().GetAttributes().SetTopologicalDimension(1);
    GetOutput()->GetInfo().GetAttributes().SetVariableDimension(1);

    vector < string > labels(4); 
    labels[0] = labelNames[0];
    labels[1] = labelNames[1];
    labels[2] = labelNames[2];
    labels[3] = labelNames[3];
    
    GetOutput()->GetInfo().GetAttributes().SetLabels(labels);
}


// ****************************************************************************
//  Method: avtTopologyFilter::ReleaseData
//
//  Purpose:
//      Release the problem size data associated with this filter.
//
//  Programmer: Akira Haddox
//  Creation:   July 1, 2003
//
//  Modifications:
//
//    Hank Childs, Fri Mar  4 08:12:25 PST 2005
//    Do not set outputs of filters to NULL, since this will prevent them
//    from re-executing correctly in DLB-mode.
//
//    Hank Childs, Fri Mar 11 07:37:05 PST 2005
//    Fix non-problem size leak introduced with last fix.
//
// ****************************************************************************

void
avtTopologyFilter::ReleaseData()
{
    avtSIMODataTreeIterator::ReleaseData();

    if (cf)
    {
        cf->SetInputData(NULL);
        vtkPolyData *p = vtkPolyData::New();
        cf->SetOutput(p);
        p->Delete();
    }
}
